Evidence of secular variation in Archean crust formation in the Eastern Indian Shield

Understanding the dominant crustal accretion model in any Archean craton is the key to understanding the dominant geodynamic process responsible for early crust formation during the Hadean (> 4.0 Ga) and Archaean (4.0–2.5 Ga). The continental crust has been proposed to have formed through either horizontal/vertical accretion related to subduction or mantle plume tectonic processes. Here, the Moho depths and average crustal Vp/Vs ratios are modelled at 16 broadband stations in the Eastern Indian Shield (EIS) through HK stacking of radial P-receiver functions (PRFs). These modelled parameters are used to test both plume and subduction models, which might have played a key role in the crustal accretion of the EIS throughout the Archean. We observe a correlation between crustal age and composition within the ellipsoidal Paleoarchean cratonic domain in the Singhbhum-Odisha-Craton (SOC), which reveals an increase in age from the younger granitoid core of the SOC (with thinning of felsic crust) to the surrounding older greenstone belts (with thickening of felsic crust). A thinner mafic crust resulting from multiple magmatic events characterizes the neighbouring Meso-Proterozoic Chotanagpur Granitic Gneissic terrain (CGGT). The Common Conversion Point (CCP) image of radial PRFs reveals northward subduction of the Paleoarchean SOC below the Meso-Proterozoic CGGT.


Seismic network, earthquake data and modelling
During 2013-2016, a local-cum-regional seismic network of 16 broadband seismic stations ( Fig. 1a; Table 1) with an average interstation spacing of 40 ± 15 km was operated by the National Geophysical Research Institute under the Council of Scientific and Industrial Research, India. The network data 13 enabled us to perform the H-K stacking 29 of radial P-receiver functions (PRFs) for estimating crustal thicknesses and average crustal Vp/ Vs ratios in the EIS. Each station was equipped with a 24-bit Reftek-130 recorder and a 120 s Reftek broadband sensor. A GPS receiver clock at each station was used for time tagging, while a recording speed of 50 samples per second was used to record continuous broadband data.
Here, three-component broadband waveforms of 150 good teleseismic events of m b ≥ 5.5 (with back azimuth between 38° and 300°, epicenters between 30°S and 90°N, and ray parameters ranging from 0.040 to 0.080 s/km) from 16 stations are used for the PRF study (Fig. 1c). First, a from − 5 to 60 s is considered from the 3-component broadband data, which is corrected from the instrumental correction. Then, all the windowed waveforms are filtered using a high pass filter with a corner frequency of 0.02 Hz before the computation of PRFs. Here, the iterative time domain deconvolution procedure of Ligorria and Ammon 30 with 200 iterations is used to compute radial and transverse PRFs for each event, using a Gaussian width factor, a = 2.5 (f < 1.25 Hz), where "f " marks the corner frequency of a low-pass filter. Thus, the wavelegth would be 4.4 km for f = 1.25 and an average crustal shear velocity of 3.5 km/s. Thus, the calculated PRFs with a = 2.5 can resolve layers with thicknesses of 2.2 km. Finally, we selected deconvolutions that reproduced more than 90% of the signal energy on the radial component (when convolved back with the vertical trace).
A total of 666 individual radial PRFs were computed using 2000 three-component waveforms from all 16 broadband stations (see Supplementary Figs. S1, S2). In this study, a minimum of 15 and a maximum of 60 individual radial PRFs at 16 different stations are used to estimate the Moho depths through HK stacking. The stacked individual radial P-RFs at 16 broadband stations (see Supplementary Figs. S2a-p, S3-11) reveal that there are at least three prominent phases corresponding to conversions from the Moho (P ms ) and crustal multiples (P p P ms and P s P ms + P p S ms phases). The Moho conversion (P ms ) and first crustal multiple (P p P ms ) show a positive conversion representing a velocity increase across the Moho, while the second crustal multiple (P s P ms + P p S ms ) shows a negative conversion representing a decrease across the Moho (see Supplementary Figs. S2-S11). The minimum arrival time of P ms is observed at CHI station, while the maximum arrival time of P ms is observed at SAL station (see Supplementary Figs. S2-S11). The back-azimuthal variations of PRF images at different stations along with theoretical arrivals of s-to-p conversion from Moho and crustal multiples are shown in Supplementary Figs. S7-S11. Radial PRF gather are also shown at different stations, showing arrivals of the Moho conversion and crustal multiples (Figs. S7-S11). Finally, the arrival times and amplitudes of the P ms , P p P ms and P s P ms + P p S ms phases from individual radial PRFs estimated for different back azimuths at sixteen broadband sites are used as inputs for the HK stacking of radial PRFs. In our analysis, we set w 1 = 0.7, w 2 = 0.2 and w 3 = 0.1 for some stations, while most of the stations gave better results for w 1 = 0.34, w 2 = 0.33, and w 3 = 0. 33 13 are also listed in Table 1.

Common conversion point (CCP) stacking of PRFs.
Here, we use the Funclab software 31 to perform CCP imaging of radial PRFs, using Dueker and Sheehan 32 's methodology to coherently stack p-to-s phase conversions for generating a 2-D image of impedance contrast at depth. The details of CCP imaging methodology have been discussed in Cladwell et al. 33 and the manual of Funclab 31 . Here, we used the 1-D IASP91 velocity model for the CCP imaging 34 . We have performed CCP imaging along five profiles (two along and three across the Himalayan collisional zone), whose locations are shown in Fig. 1a. The results of CCP stacking of radial PRFs along seven NE-SW and N-S trending profiles are shown in Fig. 5a-d and Supplementary Figs. S12a-d. These CCP images show lateral variations in the modelled Moho and LAB below the Eastern Indian Shield.
We compute PRFs using three frequency bands corresponding to Gaussian width factors, a = 1.0 (f < 0.5 Hz), and a = 2.5 (f < 1.25 Hz). Thus, the wavelegth would be 8, 4 km for f = 0.5 and 1.0, respectively, for an average crustal shear velocity of 4.0 km/s. Thus, the calculated PRFs with a = 1.0 and 2.5 can resolve layers with thickness of 4 and 1.75 km, respecively. The nominal vertical resolution at the Moho could be 4 km for the 1-D PRFs and 1.7 km for 2-D CCP stack 33 . While the horizontal resolution depends on the frensel zone. The fresnel zone width for Ps at 40 km is ~ 40 km. The station spacing for our stations varies from 10 to 20 km. So our station spacing of ~ 20 km could yield 50% overlap at 40 km depth, which is found to be suffiicient to image the Moho well 33 . At 10 km depth, the frensel zone width is 20 km, which could be image with 50% overlap for our stations with a 10 km spacing. But, our station spacing could not provide beller resolution for images at depths less than 10 km.
Following above-discussed CCP imaging technique, we generate crustal images by stacking radial PRFs. Our CCP image (Fig. 5a-d and see Supplementary Figs. S12a-d) shows a marked crustal as well as lithospheric thinning below the region covering the central Singhbhum granitoid (SG) and the north Singhbhum mobile belt (NSMB) and then both the Moho and Lithosphere deepen down toward the north below CGGT. Our modelling detects a crust-mantle boundary at 40-45 km depths beneath the region south of the SOC and thins to 35-40 km depth below the region between SG and NSMB and finally deepens to 42-48 km depths below the CGGT (Fig. 5a-d). Our modelling also detects a lithosphere-asthenosphere boundary (LAB) conversion at 160-180 km depth below the region south of the SOC (where we do not have any station) and then it thins to a depth of 130-150 km below the region between SG and NSMB and finally it deepens to 170-200 km depth below the CGGT, suggesting a northward subduction of the SOC below the CGGT (Fig. 5a-d; see Supplementary Figs. S12a-d). Note that our earlier studies using the PRF modelling 35 and joint inversion 14 of the PRFs and surface wave group velocity dispersion data have also shown a northward subduction of the SOC below the CGGT. , which is spatially correlated with the exposed Dalma volcanics (Fig. 3). Figure 4b, c show marked crustal and lithospheric thinning underlying the central granitoid of the SOC, which supports the model of vertical accretion above the Archean plume, affecting the upper mantle zone 13,20 . The zones of crustal and lithospheric thickening characterize the peripheral greenstone belts of the SOC www.nature.com/scientificreports/ (Fig. 4b, c). We observe that the mean Moho depths, lithospheric thickness, Vp/Vs value and Poisson's ratio are found to be 42 ± 5 km, 128 ± 5 km, 1.75 ± 0.12, and 0.25 ± 0.04 in the SOC, respectively, while they are found to be 40 ± 6 km, 119 ± 5, 1.83 ± 0.13, and 0.28 ± 0.05 in the EGMB, respectively (Table 1). In the EGMB, modeled crustal thicknesses and crustal Vp/Vs ratios vary from 33.6 km and 1.95 (at BHN) to 45.3 km and 1.69 (at KHU), respectively. These results get further support from the findings of our earlier studies using PRFs in the EIS 14,35 suggesting an average crustal velocity of (4.0 ± 0.08), (4.03 ± 0.06) and (3.9 ± 0.01) below the SOC, EGMB and CGGT, respectively, which, in turn, suggests a more mafic crust in the SOC and EGMB, than the relatively less mafic crust below the CGGT. In CGGT, the thinnest crust of 28.4 km thickness is modeled at RAM, which is located in the Gondwana Basin coalfields (Figs. 2, 3). The thickest crust at 46.7 km is modeled at HAZ, which is located near the northern end of the CGGT. In the central part of the Proterozoic, the CGGT-occupying HAZ and RAN stations show thicker www.nature.com/scientificreports/ felsic crust (Vp/Vs ~ 1.6-1.7). The modelled thinning of mafic crust (Vp/Vs ~ 1.8-2.0) at the LOH station in the western part of the CGGT could be due to the 1.6 Ga mantle plume activity associated with Dalma volcanism (Fig. 3a, b). From Fig. 3, another two regions of mafic crustal thinning are noticed below the EGMB and the region northeast of the study area in the CGGT, which could be attributed to the Rajmahal volcanism of ~ 117 Ma. The modeled crustal Vp/Vs ratios (Fig. 3b) reveal a marked lateral variation in the crustal composition across the EIS. The modelled average crustal Vp/Vs and Poisson's ratio estimates vary from 1.6 to 2.08 and 0.18 to 0.35, respectively. Christensen 36 suggested that minerals that exhibit Vp/Vs ratios higher than 1.8 include plagioclase, amphibole, pyroxene, and Fe-olivine. Furthermore, Christensen 36 has shown that Fe substitution for Mg in pyroxene and olivine also increases the Vp/Vs ratio; thus, basic compositions are expected to result www.nature.com/scientificreports/ www.nature.com/scientificreports/ in higher Vp/Vs ratios. Note that felsic rocks with intermediate to high silica contents could result in low Vp (< 6.7 km/s) and low Vp/Vs (< 1.78), while anorthosite rocks with high plagioclase contents have been shown to have relatively low Vp (between 6.6 and 7.1 km/s) and high Vp/Vs (> 1.85). However, mafic rocks with low silica contents could yield high Vp (> 6.7 km/s) and high Vp/Vs ratios (up to 1.86). Higher Vp/Vs values exceeding 1.86 could be possible in the presence of high-pressure fluids/melts in subduction zones 37,38 . The high Poisson's ratio values of 0.32 at CHI and 0.35 at LOH along with crustal and lithospheric thinning suggest the presence of crustal metamorphic fluids below the region west of the study area, spatially correlating with the region with exposed Dalma volcanics 25  Our modeling reveals remarkable seismic signatures within the ellipsoidal Paleoarchean cratonic domain. The central granitoid core region comprising a relatively younger SG is modelled to be characterized by a thinner felsic crust, while the peripheral regions are characterized by a thicker felsic crust, except the western peripheral side, which is characterized by the thinning of a mafic crust (Vp/Vs ~ 1.8-2.0) (Fig. 3a, b). The suggested melting of the oceanic mafic plateau model to create OMTG and IOG in the SOC at 3.52-3.38 Ga 23,24,28 is further supported by our results showing uniformly thin and felsic crust throughout the SOC region, which may dominate the initial crust formation in the Paleoarchean SOC. The primary signatures of earlier magmatism/volcanism episodes have been reported to be preserved as eruptive rocks and meta-sediments in the IOG and Simlipal basin 28 . Multiple magmatic pulses in a confined region have been suggested by the existing geochoronological data, where successive magmas signifying the granitoid batholiths follow the same conduit, resulting in pushing the older plutonic elements to the periphery of the ellipsoidal cratonic domain (involving Rayleigh-Taylor instabilities and gravitational overturns of the crust 13,39 ) to form the "keel and basin" structure of the graniteand-greenstone terranes at the surface 13 . This process has been best documented in the case of the Pilbara craton, Australia 38 . Further, the multiple magmatic pulses might have successively melted the crust and subsequently generated more felsic crust in the SOC (Fig. 3a, b). Furthermore, based on the low Cu in the Paleoproterozoic diamictites, it has been proposed that the upper continental crust became dominated by felsic rocks at the end of Archean 40,41 . Additionally, the preferential weathering of exposed mafic components might have led to the felsic Archean upper continental crust in the SOC containing lavas and intrusions, basalts, and ultramagnesian (Mgo > 18%) komatiites of 3.25 Ga 42,43 .
The neighbouring Meso-Proterozoic CGGT does not show any evidence of a "keel and basin" structure but rather provides support for the terrane accretion model that apparently worked during 1.6-0.9 Ga 15 . The mean Moho depth, lithospheric thickness, Vp/Vs value and Poisson's ratio are found to be 37 ± 8 km, 135 ± 5 km, 1.89 ± 0.21, and 0.30 ± 0.06 in the CGGT, suggesting that a thinner mafic crust characterizes the CGGT ( Table 1). The geochronological data have shown that the crustal block underlying the CGGT has been influenced by at least four magmatic-metamorphic events at > 2.5, 1.6-1.5, 1.2-1.0, and 0.9 Ga 44,45 . These events might have resulted in modeled mafic and thin crust in the CGGT. Furthermore, the presence of tholeiitic trends for basalt and calc-alkaline affinities of andesite, rhyolite and granite in the Gaya district of CGGT supports the idea of their generation in an island arc, subduction-related setting during the Neoproterozoic (1.0 Ga) 46 . Our 3-D structural plot (Fig. 4a-c) also depicts a north-dipping Moho as well as lithosphere beneath the EIS, suggesting northward subduction of the SOC below the CGGT. Furthermore, our CCP images along three NE-SW trending profiles (i.e. P1P2, Q1Q2 and R1R2 as shown in Fig. 5a-d) across the EIS detect a distinct zone-A (as shown in Fig. 5b-d) of crustal as well as lithospheric up-warping below the region between the SG and SSZ on the P1P2, Q1Q2 and R1R2 profiles, respectively, and this thinning model could be attributed to the 3.5 Ga 20 and 1.6 Ga 24 mantle plume episodes. This thinning model gets further support from the presence of ultramagnesian (Mgo > 18%) komatiites of Archean age (~ 3.25 Ga) in the Badampahar-Gorumahisani Greenstone Belt in the SOC 43 . These komatiites were crystallized from the hot magma, which were generated by moderate-high degree melting of mantle plume at different depths during the Archean 43 . Figures 5(b-d) also detect a zone-B between the SSZ and region north of it (i.e. CGGT), which is characterized by marked crustal and lithospheric thickening, supporting the idea of northward subduction of the Archean SOC below the Meso-Proterozoic CGGT that is in good agreement with the earlier proposed Proterozoic accretion model of the EIS during 1.6-0.9 Ga 15 . The CCP images along two N-S (A 1 A 2 and B 1 B 2 ) and one NW-SE (C 1 C 2 ) profiles (see supplementary Figs. S12a-d) also reveal similar images as seen in Fig. 5a-d suggesting a clear northward subduction of the SOC below the CGGT. Thus, our observations suggest that the imprints of the Paleoarchean crust formed above the mantle plume and its modification in the Proterozoic subduction process are still preserved in the EIS, which may be related to the secular cooling of the Earth's mantle.

Methods
Estimation and H-K stacking of P-radial receiver functions. The amplitude and time of P-to-S conversions (Ps) from the Moho and reverberations associated with interfaces below a recording site are generally modelled by the receiver function method (Fig. S1). The direct P-waves are generally impulsive on the vertical component seismogram at epicentral distances exceeding 30°, while Ps conversions dominate the horizontal components of ground motion (Fig. S1). The amplitude of conversion (Ps) and multiples (PpPms, PpSms, etc.) depends on the incidence angle of impinging P-waves and the size of the velocity contrasts at the interfaces. The www.nature.com/scientificreports/ P-wave incidence angle and ray parameter control the arrival times of conversions and multiples on the radial RF. All possible P-to-S conversions and multiples beneath a seismic station are shown in Fig. 2a-p, S1, S2a-p, and S3-S11.
Here, we study 666 good teleseismic earthquakes from 16 broadband stations (Fig. 1a, c). In this study, a minimum of 21 and a maximum of 35 individual radial RFs at 16 different stations are used to estimate the Moho depths through H-K stacking 29 of radial PRFs. Before stacking, a moveout correction is applied using the modified IASP91 global reference model 38 and a reference slowness of 6.4 s/° permitting summation of records from different distances. For the H-K stacking, we selected normal move-out corrected individual radial RFs (for different ranges of azimuths and epicentral distances for different stations). It is quite clear that the individual radial RFs show clear and sharp P-to-S conversions (Pms at 3.6-7.0 s after Pp (i.e., direct P arrival)) and some weak multiples (i.e., PpPms and PpSms + PsPms) from the crust-mantle boundary, suggesting that there is probably a clear Moho underlying the study region ( Fig. 2a-p, see Supplementary Figs. S2a-p, S3-S11). We note that prominent Pms arrivals characterise all individual radial RFs from all 16 broadband stations (Fig. 2a-p, see Supplementary Figs. S2a-p, S3-S11), while some weak Moho multiples are also noticed in all individual radial RFs. Finally, the crustal Vp/Vs and Moho depths were estimated at 16 broadband stations through the H-K stacking 29 of moveout corrected radial P-RFs over the available range of horizontal slowness and back azimuths ( Table 1).
The moveout times are a function of the P-wave slowness (i.e., the distance), crustal thickness, and Vp/Vs ratio. They only weakly depend on the absolute average crustal P-wave velocity 16 . In this method, the arrival times of the P ms , P p P ms and P s P ms + P p S ms phases are predicted using Eqs. (2)(3)(4), and weighting factors w 1 , w 2 , and w 3 in Eq. (1) are chosen to balance the contributions from different phases, as mentioned above. Using radial PRFs from events at different distances and stacking all the appropriately shifted traces gives a robust estimate for the thickness of the crust and (with larger uncertainty) for the Vp/Vs ratio.
After obtaining the move-out corrected radial receiver functions for different stations (using Eqs. 2-4), the average crustal Vp/Vs and Moho depth (z M ) were estimated using the approach of Zhu and Kanamori 29 . In this technique, a grid search over the Vp/Vs-z M space is performed with a view to obtain the (Vp/Vs, z M ) pair, which is the closest agreement with the observed P ms , P p P ms and P s P ms + P p S ms phases. It is observed that these phases are quite clear on the estimated radial PRFs for almost all the stations (see Supplementary Figs. S1a-b, S2a-p).
Zhu and Kanamori 29 outlined a simple method to estimate the crustal thickness and the average crustal Vp/ Vs ratio. They used a migration scheme for the direct Ps conversions and the crustal multiples for a set of receiver functions, assuming crustal homogeneity. They defined a quantity S(H,k) representing the weighted sum of the receiver function amplitudes at the calculated times of predicted arrivals of Ps, PpPms and PpSms + PsPms, which is expected to be maximum for a correct combination of H and k. This quantity can be written as: where w 1 , w 2 , and w 3 are weighting factors, which are chosen to balance the contribution from three phases viz., Ps, PpPms and PpSms + PsPms. r j (t) is the amplitude of the receiver function for the j th component, while t 1 , t 2 , and t 3 are the predicted arrival times of the Ps, PpPms and PpSms + PsPms phases.
The moveout times for the respective phases are given by the following formulas: where h denotes the crustal thickness, and a and b are defined as: and p is the horizontal ray parameter. Furthermore, we model Poisson's ratio from the estimated Vp/Vs values using the following relation:

Data availability
Most of the data that support the findings of this study are available in the supplementary material. (1) S(H, k) = w 1 r j (t 1 ) + w 2 r j (t 2 )−w 3 r j (t 3 )